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PROGRAMS FOR COMPUTING EQUILIBRIUM THERMODYNAMIC 


PROPERTIES OF GASES 
By Harry E. Bailey 
Ames Research Center 


SUMMARY 


This report presents a series of subroutines, written in FORTRAN IV, for 
calculating the equilibrium thermodynamic properties of any mixture of react- 
ing diatomic gases if the temperature, density, and molar concentrations of 
the species are given. One triatomic gas, carbon dioxide, may also be 
included in the mixture. Two additional subroutines permit the calculation 
of the molar concentrations of the species if the gas mixture is air. A 
sample main program is included which uses all these subroutines to compute 
the equilibrium composition and thermodynamic properties of air for a given 
density over a range of temperatures. 


INTRODUCTION 


Solutions to problems involving the flow of a chemically reacting gas 
require the computation of the equilibrium thermodynamic properties of the 
gas mixtures over wide ranges of temperature and density. These properties 
are available in tabular form for air (refs. 1, 2, and 3), in graphical form 
for carbon dioxide (ref. 4), and for three nitrogen carbon dioxide mixtures 
(ref. 5)* However, neither the tabular nor graphical form is convenient for 
use in electronic machine calculations without elaborate and tedious curve 
fitting of the data. Therefore, the subroutines presented here should prove 
useful for the solution of many problems in the flow of a chemically reacting 
gas. 


Five different subroutines are presented in this report. In addition, a 
main program designed to compute the equilibrium thermodynamic properties of 
air is presented to illustrate the use of the five subroutines. Two of the 
subroutines, SPECIE and ENTH, can be used to obtain the pressure, enthalpy, 
specific heat at constant pressure (with frozen composition), entropy, gas 
constant, and equilibrium constants of the components of a gas mixture. The 
various chemical and physical constants needed as input by subroutine SPECIE 
may be obtained from references 6 and 7 . A third subroutine, RATCON, gives 
the reaction rates for eight chemical reactions in air. This subroutine 
would require modification if some other gas mixture were used. Two addi- 
tional subroutines, GUESS and ITERA, calculate the equilibrium molar concen- 
trations of each species for air at a given temperature and density. 

Although these subroutines would require major modification for any mixture 
other than air, they do provide a suitable basis from which similar programs 
for other mixtures could be written easily. 


I 


Each individual species is assumed to behave as an ideal gas. The 
thermodynamic properties of the diatomic species are approximated by the 
harmonic-oscillator rigid-rotator model with rotational and vibrational 
constants appropriate to the lowest electronic state. Computational results 
are given for two air models. Model A contains the species N 2 , 0 2 , N, 0, NO, 
N0 + , 0 + , F 4 ", and e“. Model B contains the same nine species as model A plus 
two additional species, 0 ++ and N ++ . 


SYMBOLS 

c specific heat at constant pressure (with composition frozen) 

Cp^ specific heat at constant pressure for the ith species 

energy of the nth energy level 

free energy of the ith species 

g^ degeneracy of the nth energy level 

H total enthalpy of the mixture 

EL enthalpy of the ith species 

enthalpy of formation of the ith species 
u i 

KL equilibrium constant for the ith reaction 
PL molecular weight of the ith species 

n^ number of atoms in the ith species 
P total pressure of the mixture 

P q standard pressure of one atmosphere 

p^ partial pressure of the ith species 

R universal gas constant 

S total entropy of the mixture 

entropy of the ith species 
T temperature, °K 

T characteristic rotational temperature 

characteristic vibrational temperature 
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ZN mass concentration of nitrogen atoms in cold mixture 
ZO mass concentration of oxygen atoms in cold mixture 
3. sum of 3. . values over the index 
3. . difference of the stoichiometric coefficients of the jth species in 

ij i -i • • -> i * 




the ith reaction 
concentration of the ith species in moles/gm 


SAMPLE MAIN PROGRAM 


This is a typical program (see appendix A) which uses all the following 
subroutines to compute the thermodynamic properties of equilibrium air. It 
is intended as a guide for the use of the subroutines SPECIES, ENTH, RATCON, 
GUESS, and ITERA. It by no means exhausts the possibilities for their use. 

The sample program requires an input AT (DELTA), a density (RHO), and 
nine species concentrations for the cold mixture of nitrogen and oxygen (g(n)). 
The program will then call the appropriate subroutines for the computation 
of the thermodynamic properties pressure, gas constant, enthalpy, entropy, 
specific heat at constant pressure, equilibrium constants, reaction rate con- 
stants (as defined under subroutine RATCON), and the species concentrations. 
These computations will be performed at the given density for 25 temperatures 
starting at AT and ending at 25 AT. A sample output is given for air at a 
density of 0.01288 gm/cm 3 and at 4 temperatures from 1000° K to 4000° K. 


SUBROUTINE SPECIE 


This subroutine (see appendix B) will compute the pressure (P), 1 
enthalpy (h( 21)), specific heat at constant pressure (with composition 
frozen) (CP(21)), entropy (FE(21)), equilibrium constants (BK(j)), reaction 
rate constants (AK(j)), and the gas constant (GCONST) for a given temperature 
(T) , density (R), and the molar concentrations (GA(j)) of each species. On 
the first entry to the subroutine, all necessary chemical and physical 
constants are read from cards and stored in the program. Provision has been 
made for the computation of 20 species and 40 reactions. 

This subroutine uses the following equations to compute the thermodynamic 
properties of the mixture from the thermodynamic properties of the components. 


^■Characters in parentheses are those used in the FORTRAN listing in the 
appendixes. There is an unavoidable inconsistency between these characters 
and those defined under Symbols. 
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The enthalpy is given by 


\ 


_H. _ 
RT 


H i 

7 i EF 


( 1 ) 


The specific heat at constant pressure (with frozen composition) is 
given by 


■ I 


P /_> 'iTi 

The entropy of the mixture is given by 


1= In 


Hi F ± pt 
RT RT 


The total pressure of the mixture is 

P = pRT 


I 7 i 


The equilibrium constant for each reaction is 

ln K i = "X P i«3 " P i ln ( 82 ' °5576l T) 

The partial pressure of each species is computed from 

P. = P 

l £7 . 

' i 

The entropy of each species is computed, in this subroutine as 


( 2 ) 


( 3 ) 


( 4 ) 


( 5 ) 


( 6 ) 


s. 

1 




( 7 ) 


The data that must be input to this subroutine are the number of 
reactions (IL), the number of species (JL), the stoichiometric coefficients 
for each reaction (B(l,J)), the number of electronic energy levels for each 
species (NL(j)), the energy (E(j,N)) and degeneracy (G(<J,N) ) of each of these 
electronic energy levels, the number of atoms in each molecule (ENJ(J)), the 
vibrational (TV(j)) and rotational (TR) temperatures of each diatomic species, 
the heat of formation (ENTA(j)) of each species, the molecular weight (EMWT(J)) 
of each species, and finally the percent by volume (GX(j)) of each species in 
the cold mixture. All this information is stored in the subroutine and used 
on all subsequent entries for the computation of the thermodynamic properties. 

Each time the subroutine is called it must be given the temperature (T), 
in degrees Kelvin, the density (R) in gm/cm 3 , and the mass concentrations 
(GA(l)), in moles/gm. The subroutine will then compute all the quantities 
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listed in the first paragraph of this section. In addition, the enthalpy 
(h(J)), entropy (FE(j)), and specific heat at constant pressure (CP(J) ) for 
each individual species (J = 1, JL) will be computed. 

Subroutine SPECIE uses two subroutines, ENTH and RATCON, which are 
discussed in subsequent sections. 

Sample input data required by subroutine SPECIE are shown in appendix F. 
The input data are shown for model A and model B. 


SUBROUTINE ENTH 


This subroutine (see appendix C) will compute the enthalpy (H), specific 
heat at constant pressure (CP), and free energy (FE) of a given species for 
a given temperature (lX). For each species, this subroutine must be given 
the heat of formation (ENTAX) , the number of atoms in one molecule of the 
species (ENJX), the vibrational temperature (TVX), the number of electronic 
states (NLX) , as well as the energy (EX) and degeneracy (GX) of each elec- 
tronic state and a constant (EMWTX) which depends on the molecular weight 
of the species and the rotational temperature of the species. 


The computations performed by this subroutine are based on the harmonic- 
oscillator rigid -rotator model for all diatomic molecules. It is also pos- 
sible to compute the thermodynamic functions for one triatomic molecule, 
carbon dioxide. In this computation it is assumed that there is only a 
ground electronic state for the carbon dioxide molecules. 


The equations used by this subroutine for the computation of the 
thermodynamic functions of each individual species are as follows. The free 
energy is computed from 


F i 

RT " ~ ln 


e "E n /T _ [( n _ i) + 2.5 ]Zn T - 1.5 In M. + 3-6649516 

°n v l i 


(n ± - 1) In (l - e" Ty / T ^) + (n ± - l)ln T r + 


The enthalpy of each species is 


H i _ r, 

m~ ^ n i 


1) + 2 5] + i. Sg n E n e ~ En//T + ( n i ~ ^ T v + ^^£i 

1} 2 ’ 5J T “ -En/T 3 Tv/T 3 RT 


2gne' En/r T(e Tv/T - 1) 


(8) 


( 9 ) 
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The specific heat at constant pressure (with frozen composition) is 


= [(n. - 1) + 2.5] + (n, - 1) fe 


a e T v /T 


(e T vA . x 

! (^n E / e ' En/T ) - kV^ 


T 2 


2g n e 


-E n /T 


( 10 ) 


SUBROUTINE RATOON 


This subroutine (see appendix D) will compute the rate constants for a 
given temperature according to the formulas and constants of references 8 and 
9 for the reactions listed below. 


>2 + O 2 

-► 

20 + o 2 

k f 

(11) 

Ig + N 2 


2N + N 2 

K f 

(12) 

NO + M 


N + 0 + M 

% 

(13) 

N + 0 2 

3± 

NO + 0 


(14) 

0 + n 2 

4- 

NO + N 

k f 

(15) 

0+0 

* 

4* - 

0 + e 

k b 

( 16 ) 

N + N 


N + + e" 


(IT) 

N + 0 

— ► 
4- 

N0 + + e" 

% 

(18) 


The symbol Kjp after a reaction indicates the forward rate constant; 
Kg indicates the backward rate constant. 

This routine requires only the temperature as input. It is a simple 
matter to alter this subroutine to give the rate constants for any reaction 
desired. 


SUBROUTINES GUESS AND ITERA 


For any mixture of nitrogen and oxygen^ two additional subroutines 
(see appendix E) are available that permit the computation of the equilibrium 
species concentrations at a given temperature and density. These two 
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subroutines require as input the temperature (T), the density (RHO), the 
atom concentrations of the cold mixture (ZN,Z0) , and eight equilibrium 
constants (EK(l)). 

The first subroutine GUESS estimates the relative magnitude of the 
various species concentrations. The second subroutine ITEM then refines 
these guessed values until they satisfy the three conservation equations, 
that is, conservation of oxygen atoms, nitrogen atoms, and charge, as well as 
six equilibrium equations. These two subroutines require a block of common 
storage called /C0M4/ which contains arrays for 20 equilibrium constants 
(EK(20)), 20 mass concentrations (B(20)), 20 first guesses (c(20)) of the mass 
concentrations, the density (RHO), and the cold mass concentrations of 
nitrogen atoms (ZN) and oxygen atoms (ZO) . 

The equations that must be solved by subroutine ITEM are listed below. 


20 = 7 o + 27 o a + 7 N0 + 

7 N0 + + 7 0 + 

(19) 

2N = 7 n + 27 n 2 + 7 N0 + 

7 N0 + + 7 F f 

(20) 

e- = 7 N0 + + 7 0 + + 7 N+ 


(21) 

P7 q 2 

Ki - 0 


(22) 

7 0 2 



k 2 = p7 * 2 


(23) 

0J 



Ks - 


(24) 

7 NO 



Ke - 


(25) 

’’O 



k 7 = 


(26) 

7 n 



K S . h£»oi 

7„7 m 


(27) 


The convergence of the iteration process used in subroutine ITEM 
requires that equation (19) be solved for the largest of /q, 7 q 2 , or 7 q+, and 
that equation (20) be solved for the largest of 7^, y^ z , or 7™+. The remain- 
ing 7^ values are then found from equations (22) through (27). The value 
of 7 g _ is always found from equation (21). The path taken in subroutine 
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ITERA depends therefore on the relative magnitudes of jq , yq 2 , 79+ and of 
7 jj, 7n 2 > an< ^ ^N + * '^ ie function of subroutine GUESS is to provide an initial 
guess of these relative magnitudes. This initial guess is based on the 
temperature and density of the mixture and on the relative values of the 
equilibrium constants Ki and K s for oxygen and K 2 and K 7 for nitrogen. 

Based on these guessed values, subroutine ITERA then computes the 
species concentrations by iterating equations (19) through (27) until all the 
equations are satisfied. At each step in the iteration the newest value of 
each concentration is computed as the average between the value obtained on 
the previous iteration and the value computed on the present iteration. The 
iteration is terminated when the change in the values for two successive 
iterations is less than a specified amount (10 -7 for model A and 10 -5 for 
model B) and when all the equations are satisfied to within a specified 
amount (5xl0 _s for model A and 10" 6 for model B). 

Appendix E gives two listings for subroutine ITERA and two listings for 
subroutine GUESS. In each case one listing is for model A (contains no 0 ++ 
or N ++ ) and one listing is for model B (contains 0 ++ and N ++ ). 


POSSIBLE USES OF THE PROGRAMS 


A series of programs for computing the thermodynamic properties of a 
variety of gases is presented. In addition, two of the programs permit the 
computation of the equilibrium species concentrations of air or any other 
mixture of nitrogen and oxygen, as well as the thermodynamic properties of 
the mixture. 

The programs are presented in a form which can be modified easily to 
suit the user's needs. For instance, the thermodynamic properties of any 
diatomic or monatomic gas may be obtained simply by altering the input data 
to subroutine SPECIE. Subroutine RATCON can be easily altered to give any 
rate constant desired. Subroutines GUESS and ITERA may be rewritten to give 
the equilibrium concentrations of gas mixtures other than air. All the sub- 
routines as presented may be used in a new main program to compute the 
equilibrium flow properties behind a normal shock wave. 


DISCUSSION OF RESULTS 


The main program described in this report has been used to compute the 
thermodynamic properties of air. Two different models have been used. In 
model A the species assumed to be present are Ng, 0 2 , NO, N, 0, N0 + , N + , 0 + , 
and e“. Model B contains N ++ and 0 ++ in addition to the species in model A. 
The values obtained are compared to the values presented in references 2 and 
3 (table I) and to the values presented in reference 1 (table II). 
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The dimensionless enthalpy, H/RT, computed for model A agrees with the 
, values presented in reference 2 to within 1 percent except at 15,000° K, where 
the error rises to 1.39 percent at 100 times standard density and to 8.46 per- 
cent at 10 -7 times standard density. The dimensionless entropy, S/R, computed 
for model A, agrees with that computed in reference 2 to within 0.5 percent 
except at 15,000° K and 10 -7 times standard density where the error is 
2.54 percent. The pressure (in atmospheres), for model A, agrees with that 
computed in reference 2 to within 3 percent except at 15,000° K and 10 “ 7 times 
standard density where the error is 3*27 percent. 

The increase in error at 15,000° K for enthalpy, entropy, and pressure is 
attributed to inadequacies in model A. The neglect of species IT*' 4 ' and 0 ++ 
begins to be felt at 15,000° K. That this is indeed the case may be seen by 
comparing the results computed for model B to the results of reference 2. At 
15 ,000° K and p/p Q = 10 “7, the error in the dimensionless enthalpy drops from 
8.46 to 0.13 percent, the error in dimensionless entropy drops from 2.54 to 
0.25 percent, and the error in pressure drops from 3*27 "to 0.28 percent. 

The highest temperature available in references 2 and 3 is 15*000° K. 
Therefore, for comparisons at higher temperatures, it is necessary to use the 
results of reference 1. The comparison between values in reference 1 and the 
present calculations is presented in table II. Again, comparison is made with 
both the model A and model B air of the present report. At 18,000° K and 
25,000° K and at p/p Q = 10“ 6 , the errors in model A are seen to be enormous. 
However, for model B, which includes N ++ and 0 ++ , the errors at these 
points are drastically reduced. In fact, they become comparable to the errors 
in the low temperature region. 

CONCLUSIONS 


A series of subroutines written in FORTRAN IV is presented. These 
programs may be used to compute the thermodynamic properties at a given den- 
sity and temperature of any diatomic gas plus the one triatomic gas, C0 2 , if 
the molar concentrations are given. In the special case of air, additional 
subroutines are included that permit the computation of the molar concentra- 
tions themselves. The thermodynamic properties of two models for air are 
compared to the results obtained by other investigators. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., Dec. 28, 1966 
129-01-08-20 
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APFEM)IX a 


$I8FTC MAIN NCREF 
CHAIN CHECKOUT BAILEY 

C0MMCN/CCM4/EK<20> ,8(20) «C( 20) ,RHO,ZN,ZO 

01 HE NS I ON AK(20)«BK(20),H(21),CP(21),S(21),G(21) 

RE AO I 5 » 1 2 ) NHA X 

12 FCRHAT ( I 3 ) 

RE AD ( 5 , 2 ) DELTA 

2 FCRHAT ( F 15 • 6 ) 

REAC 15,4 ) (GIN), Ml, 9) 

4 F0RMATI5F15.6) 

ZN-2 • *G ( 5 ) 

ZC = 2 • *G I 4 ) 

I RE AC I 5 , 2 ) RFC 
TEMP*0 • 

WRITE<6,13) 

13 FORMAT! 12H1TEMPERATURE/79H CENS l TY , PRESSURE , GAS CONSTANT ,ENTHA LPY , 
1ENTR0PY, CONSTANT PRESSURE SPECIFIC HEAT/9H ENTHALPY/ 8H ENTROPY /3 2H 
2 CONSTANT PRESSURE SPECIFIC HE AT/2 1H EQUILIBRIUM. CONST AN T/ 1 AH REAC 
3 T I ON RATE/ 19H HASS CCNCENTRATI0N/34H 0, N, E- , 02 , N2 , NO ,N0^ ,0 ♦ , N«- ,0* + 
4 ,N ♦ ♦ ) 

DO 3 M*1 , NH A X 
TEHP*TEHP+DELTA 

CALL SPECIE! TEMP,RHO.PRESS,G*H, CP, S , AK , BK, GC0NS1 ) 

DO 0 N* 1 , l i 
IF (BK(N)-82. ) 9,9,10 

9 EK(N)*10.*«BKIN) 

GO TO 8 

10 EK I N I * 1 • E + 36 

8 CONTINUE 

CALL GUESS 
CALL ITERA 
DO 11 N*i,ll 

II GIN)*B(N) 

CALL SPECIE!TEHP,RHO,PRESS,G,H,CP, S,AK,BK,GCONST) 

WRITE <6, 5) TEMP , RHO, PRESS, GCONST, HI 211 • SI 21), CPI 21 ) 

5 FORMAT! 1H0,E15.6/10X,6E15. 6) 

WR I T E 1 6 , 6 ) (H(N),N*l,ll» 

WR I T E ( 6 , 6 ) (SIN), N* 1,11) 

WR I T E ( 6 , 6 ) I CP IN), N* 1,11) 

WRITEI6,6) ( BK ( N ) , N= 1 , 10 ) 

WRITEI6.6) ( AK ( N ) , N 3 1 , 10 ) 

WR I TE ( 6 , 6 ) IGIN) ,N*1 , 11 ) 
t FCRHAT ( IH ,11EII.A) 

3 CONTINUE 
WRITER, 7) 

7 FORMAT! 1H0, 10HENC OF RUN) 

GC TC 1 
END 
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APPENDIX C 


SIBFTC TC 1210 NGREF 

SUBROUTINE ENTH(TX,ENTAX,ENJX*TVX,EMWTX,EX»GX*NLX*H,CP»FE) 
DIMENSION EX < 20 ) *GX l 20 ) 

NX-IFIX(ENJX) 

GO T0(1, 2,3), NX 

1 SUM1-0. 

SUM2*0 • 

SUM3*0. 

DO A N=1 *NLX 
XX=GX(N)*EXP(-EX<N)/TX> 

SUM l*SUNl+XX 
SUM2*SUM2+XX*EX(N) 

A SUM3=SUM3*EXIN)*EX(N)*XX 

FE*-AL0G($UM1>-2.5*AL0G(TX)-EMWTX+ENTAX/TX 
H=2.5+SUM2/(TX*$UM1)+ENTAX/TX 
CP=2.5+(SUM1*SUM3-SUM2«*2)/ l l TX«SUM 1 ) **2 ) 

GO TC 5 

2 SUM 1*0 • 

SUM2*0* 

SUN 3*0. 

DO 6 N= 1 » NLX 
XX=GX(N)*EXPI-EX(N)/TX) 

SUMl^SUMUXX 

SUH2 s SUM2+XX*EX!N) 

6 $UM3*SUM3+EX(N)«EX(N)*XX 

FE*-AL0G(SUM1)-3.5*AL0G(TX)-EMWTX»ENTAX/TX 
H-3.5+$UM2/(TX*SUMl l+ENTAX/TX 
CP=3.54(SUM1*SUM3-SUM2**2)/(CTX«SUM1 >**2> 

IF ( TVX/TX-80 *0 ) 7,7,8 

8 V I BH*0 • 

VI BCP*0 • 

R2-0. 

GO TO 9 

7 A1 =EXP I TVX/ T X ) 

R2*AL0GU.-EXP(-TVX/TX ) I 
A11=AI-1. 

VIBH*(TVX/TX)/A11 

VIBCP*(TVX/TX)**2 

VIBCP*tVl8CP/All)*Ul/All) 

9 CP = CP* V I BCP 
H=H+VI BH 
FE-FE+R2 

GO TO 5 

3 VIB1-1932.1/TX 
V I 62*960 • 1/TX 
VIB3*3380./TX 
VIB11*1./(1,-EXP(— VIB1) ) 

VIR22~l*/( 1*-EXP(— VIB2) ) 

VIB33~l*/I 1*-EXPI— VIB3) ) 

VIB1C* VIB1*«2*EXPI VIBl)/< EXP t VIED- 1.)*«2 
VIB2C S VIB2**2*EXP(VIB2)/(EXP(VI62)-1«)*«2 
VIB3t=VI83**2*EXPIVIB3)/<EXP(VI83)-l.)**2 
VIB1*VIB1/(EXP(VIB1)-1.) 

VIB2*VIB2/(EXP(VIB2)-1. ) 

VIB3=VIB3/(EXP« VIB3I-1. ) 

H=3,5 + VIBl*2.*VIB2 + VIB3-fENTAX/TX 

FE = -3.5*ALCGtTX)-At0GIVIBU*VIB22*VIB22»VIB33)-EMWTX+ENTAX/TX 
CP“3*5+VIB1C*2«*VIB2C+VIB3C 
5 RETURN 
END 


APPENDIX D 


$ I BFTC TC 1 208 NOREF 

SUBROUTINE R ATCONI AK , TX ) 

DIMENSION AK(IO) 

T1 5*TX*SCR T < TX ) 

AKU)=3.6E*21«EXP(-59380./TXI/T15 

AK(2)=l.5E+20/T15 

AK(3)«5. 18E*2l*EXP( -75490. /TXI/T15 
AKU)*l.GE+12*EXPC-3120./TX)*SCRTtTX) 
AKI5)*5.0E*13*EXPI-36016./TX> 
AKI6)-6.0E*24/(T15»TX) 

AK ( 7 ) =AK 1 6 ) 

AK(8)-U8E + 21/TX5 

RETURN 

END 
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ENC 


MODEL A 


i 


tlBFTC TCI 206 NCREF 

SUBROUTINE ITERA 

C0HNCN/C0P4/EK120) ,BI20 ),C(20) ,RHO,ZN,ZO 
DATA DELl,DEL2/l.E-07,5.E-08/ 

DO 9 N=l,200 
IF(B14)-BI1>) 1,1,2 

2 CM )-( ZC-B ( 1 1-6 (6)-B (71-8(8)1/2. 

BM )*l 8 (4 ) *C 14 H/2. 

Cl 1)*SQRT(B(4)*EK( 1 J/RHQ) 
Bll)>l8ll)K(l)i/2. 

CI8MIBI 11/8(3) )«(EKt6l/RHO) 
B(6)»(B(0)+C(8) )/2. 

NTEST1-1 
GO TO 3 

1 IF(B(1)-B(8M 4,4,5 

5 Cll WC-2.»B(4)-ei6)-Bl7)-B(e> 
Bdl-IBClUCm )/2. 

C(4)*RH0*B{ l)*B(l)/EK( 1) 

B<4 )*IB14)+C 14 ) )/2. 

C(8)*(B( U/813) )*(EK(6)/RHO) 

8 ( 8 ) * I B ( 8 ) +C ( 8 ) )/2« 

NT EST1*2 
GO TO 3 

4 Cl8)*ZQ-B(l)-2.*B(4)-8(6)-Bl7) 

B18)*IB18>*C(8) )/2. 

Cl 1 )-RHC«B(3)«Bte)/EKl6) 

B ( 1 > *1 B f 1 MCll) )/ 2 . 

C(4)*RH0*8( 1 ) *8 l 1 ) / EK ( 1) 

B14 »-(Bl4)+Cl4 ) )/2. 

NTEST1-3 

3 IF ( B 1 5 ) -B ( 2 )) 6,6,7 

7 Cl5)«lZN-8l2)-B(6)-Bm-Bl9))/2. 

B( 5 ) *(B l 5) +C 15 I )/2. 

C(2I*SQRT(B(5J*EK(2 I/RHOI 
B(2)*(B(2)+CI2) )/2. 

C19)-IB12I/B13) ) «« EK 17 )/RHO ) 
B(9)*IB19)+C19) J/2* 

NTEST2-1 
GO TO 8 

6 IF f B 12 )-BI 9 ) ) 31,31,32 

32 Cl2)-ZN-2.*ei5)-8l6)-B(7»-B(9) 

B(2)*(B(2)+C(2) ) /2 • 

Cl 5 I *RHO*B (2)«B(2)/EK(2) 
B(5)*lB(5)+Cl5))/2. 

Cl 9 ) *(B 1 2 ) / G ( 3 ) J •! EK 17 I /RHO ) 
B(9)»IBI9IKI9)|/2. 


NTEST2*2 
GO TO B 

31 C19)*ZN-B(2)-2.*B15I-B(6)-B17J 

B(9)*(B(9)+C(9) )/2. 
C12)*RH0»813)*B19)/EK17) 
B(2)*(B(2)+C(2))/2. 

C< 5)*RHC«B(2 )*B12)/EK12) 

B(5)»(B(5I+C(5J 1/2. 

NT E ST2 3 3 

8 C(3)*B(7)+B(8)+B(9) 

B(3)*(B(3MC(3) )/2. 
C16)»RH0*BU)*B12)/EK(3) 

BI6)*IBI6)+C(6) )/2. 

Cl 7)*B(1)*B12) •EK18)/B(3) 
B(7)*(Bm+Cl7U/2. 

10 IFIABS I (Cm-Bll) l/eiil )-DELl ) 11.11,9 

11 IFIABS l ICl 2J-B12I )/Bl2) )-DELl ) 12,12,9 

12 IFIABS (ICC3)-B(3I I/BC3I ) — DELI ) 13,13,9 

13 IFIABS l IC14I-BI41 )/B(4) i-DELl ) 14,14,9 

14 IFIABS l (C15)-B(5) )/B(5) )-DEL 1 ) 15,15,9 

15 IFIABS l CC(6I~B(6) 1/8(61 l-OELl ) 16,16,9 

16 IFIABS l f C I 7 1 — B 171 ) / B ( 7 I J-CELl ) 17,17,9 

17 IFIABS UCm-Bie) )/B(8) )-DELl ) 18,18,9 

18 IFIABS l (C19J-B19) )/B(9) ) - DEL 1 ) 19,19,9 

19 DWC-Bl i)-2.«B(4)-B(6)-B(7 1-0 1 8) 
02*ZN-B(2)-2.*6<5)-B(6)-B(7)-BI9) 
03*BI3I-B(7)-B(8)-B19) 

04*B(4I*EKl 1 )-RHC* B ( 1 ) *B ( 1 ) 

05*8(5) «EK12 )-RHC«Bl 2)*Bl 2 ) 

06*B<6) *EK( 3 )-RHC*8( l)*Bl 2 > 

07* B( 1 ) •EKl 6 )-RhC*Bl 3 ) *B ( 8 ) 

08 *B l 2 ) *EK m-RHOBO) *B19 ) 

D9*B(1) *B(2 )*EK(8)-B(3)»B17) 

IFIABS (0U-CEL2 ) 22,22,9 

22 IFIABS ( 02 ) -DEL2 ) 23,23,9 

23 IFIABS (03 ) -0EL2 ) 24,24,9 

24 IFIABS l 04 ) -CEL2 ) 25,25,9 

25 IFIABS ( 05 ) -CEL2 ) 26,26,9 

26 IFIABS ID6J-CEL2 ) 27,27,9 

27 IFIABS l D7 ) -DEL2 ) 28,28,9 

28 IFIABS I D8 ) -CEL2 ) 29,29,9 

29 IFIABS (09J-0EL2 ) 30,30,9 

30 GO TC 33 

9 CONTINUE 

MRITE(6,34)NTEST1,NTEST2, (Bl J) ,J*1.9) 

34 FORMAT! 1H0.22HN0 CONVERGENCE NTEST1-13, 
7HNTEST2-I3/IH ,1P9E14.6) 

CALL EXIT 
33 RETURN 
END 
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Letters in ( ) indicate gas model. 
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